Network pruning.

Contents.

Introduction.

We start from the Gene Regulatory Network (GRN) detailed in the figure below. In this notebook we show:

initial_network.png

Diagram of the stem cell developmental network. It is made up of 52 nodes, which represent different genes. The connections may represent either activations (arrows), or repressions (perpendicular, T-shaped).

Imports.

Each package can be installed with the commented lines above every using instance.

Parameters.

We state the parameters used to analyze the system of ODEs that represents the GRN:

$\dot{x}_i = F_i,\;\; i=1,2,...,52, \label{eq:52}$,

where $F_i$ is given by

$F_i = \sum_{j\in A_{ij}}{\frac{ax_j^n}{S^n + x_j^n}} + \sum_{j\in B_{ij}}{\frac{bS^n}{S^n + x_j^n}} - kx_i$,

where $S$ controls the strength of the regulatory interactions by allowing to tune the threshold value above on which the activator/repressor has an effect. The Hill coefficient $n$ determines the sigmoidal function's steepness. The self-degradation constant is $k$. The repression constant is $b$. The activation constant is $a$. In the summation subindex, $A_{i,j}$ represents the set of activator nodes acting on $i$, while $B_{i,j}$ represents the set of repressor nodes acting on $i$. The first term in the equation represents self-degradation, the second term the activation from node $i$ to node $j$, and the third term the repression from node $i$ to node $j$.

m52 is a list of lists that describes the interactions between nodes in the network.

For instance, the first 10 interactions in the network are shown below. The first element of each list represents the node origin of the interaction, and the second represents the node target. The third element is either 1 for activation or -1 for repression.

Besides, we also need to declare some parameters to repeatedly solve the system of ODEs:

It is important to explore the whole space of initial conditions. The minimum possible concentration is 0. We calculate the maximum concentration $x_i$ in the function $F_i=0$ by setting $x_j \rightarrow +\infty,\; \forall j\neq i$ and $x_j \rightarrow 0,\; \forall j\neq i$. The former case yields

$x_i = (\sum{a})/k$,

while the latter yields

$x_i = (\sum{b})/k$,

thus, the maximum possible concentration of $x_i$ corresponds to the addition:

$x_i = (\sum_{j\in A_{i,j}}{a} + \sum_{j\in B_{i,j}}{b})/k$,

Our initial conditions will be random concentrations, uniformly distributed between 0 and 2.88---one concentration per node.

Data and functions.

The following function uses the list m52 (or sparse matrix) to define the system of ODEs.

The following functions are used to:

The following variable names_ordered includes the names of each of the 52 genes. In the matrix m52, node 1 is OCT4, 2 is SOX2, 3 is NANOG, etc. GATA6 is node 16.

Function for plotting:

Bistability of the original GRN.

With our functions it is very simple to repeatedly solve the system of ODEs a number of times:

Once we obtain the solution, we can calculate the attractors of a node. The outputs are:

Besides, we can also easily plot the solutions of a node:

Pruning the original network.

We'll begin pruning by selecting the downstream nodes (those that receive connections but do not send them). Then, from the resulting network, keep recursively selecting the downstream nodes until none are left.

Prune the network by checking whether each connection in m52 corresponds to a downstream connection. If so, delete the connection in m52_prunned, which is initially a copy of m52.

Now we repeat the process iteratively, until the network cannot be reduced anymore.

The code above does note identify as downstream nodes that have a self loop.

We check if there are any:

Once the downstream nodes that have a self-loop have been removed there may be other newly downstream nodes.

We now have a new sparse matrix to describe the system of ODEs of the pruned network. The problem is that we have kept the same numeration of nodes as in the pre-pruned network. We need to re-write the order of the nodes.

We call m26 to the resulting sparse matrix (because it has 26 nodes, instead of the 52 of the original):

One may save the resulting list of lists m26 to a text file by copy-pasting the list Julia prints:

Get the pruned network in Cytoscape format:

Uncommenting the lines below, one can import the network obtained here into a format suitable to visualize in Cytoscape.

The reader needs only to uncomment the cells below, and, in the functions open(), add a directory to save a CSV file.

The resulting network looks like this after importing to Cytoscape:

network_w_downstream_reviewed.png

Each node is colored depending on whether it is upstream (i.e., source of connections butnot target, in green), downstream (i.e., target of connections but not source, in red) or neither (source and target of connections, in yellow).

Ignoring the downstream nodes, we obtain:

network_pruned_reviewed.png

Bistability of the 26-node network.

New function to describe the system of 26 ODEs:

Besides, we can also easily plot the solutions of a node:

Knock-out genes experiment.

We reasoned that the 26-node network may be further reduced. Hence, we now perform knock-out experiments on each gene in the pruned network, and keep save a list of those that have no influence on the bistability (within a margin error eps).

The nodes that are relevant for keeping bistability are:

We can save and print their names:

We hypothesize that removing all of the nodes (at the same time, nodes at the left of figure below, in cyan) that have no influence on bistability keeps the original bistability unchanged. We can check this by simulating the system of 13 ODEs and comparing to the original solutions.

network_w_downstream_blue.png

First, we build the sparse matrix that describes the system of ODEs in this 13-node network:

We call m13 to the resulting sparse matrix (because it has 13 nodes, instead of the 52 of the original):

One may save the resulting list of lists m13 to a text file by copy-pasting the list Julia prints:

The resulting network is:

13_nodes_clearer.png

Bistability of the 13-node network.

Now that we have the sparse matrix m13 we can solve the corresponding system of ODEs.

Besides, we can also easily plot the solutions of a node: